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Abstract 

Proper  orthogonal  decomposition  (which  is  also  known  as  the  Karhunen 
Loeve  decomposition)  is  a  reduction  method  that  is  used  to  obtain  low 
dimensional  dynamic  models  of  distributed  parameter  systems.  Roughly 
speaking,  proper  orthogonal  decomposition  (POD)  is  an  optimal  tech¬ 
nique  of  finding  a  basis  which  spans  an  ensemble  of  data,  collected  from 
an  experiment  or  a  numerical  simulation  of  a  dynamical  system,  in  the 
sense  that  when  these  basis  functions  are  used  in  a  Galerkin  procedure 
will  yield  a  finite  dimensional  system  with  the  smallest  possible  degrees 
of  freedom.  Thus  the  technique  is  well  suited  to  treat  optimal  control 
and  parameter  estimation  of  distributed  parameter  systems.  In  this  pa¬ 
per,  the  method  is  applied  to  analyze  the  complex  flow  phenomenon  in 
a  horizontal  chemical  vapor  deposition  (CVD)  reactor.  In  particular,  we 
show  that  POD  can  be  used  to  efficiently  approximate  solutions  to  the 
compressible  viscous  flows  coupled  with  the  energy  and  the  species  equa¬ 
tions.  In  addition,  we  also  examined  the  feasibility  and  efficiency  of  POD 
method  in  the  optimal  control  of  the  source  vapors  to  obtain  the  most 
uniform  deposition  profile  at  the  maximum  growth  rate.  Finally,  issues 
concerning  the  implementation  of  the  method  and  numerical  calculations 
are  discussed. 
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1  Introduction 

Chemical  vapor  deposition  (CVD)  processes  use  a  chemical  reaction  in  the  gas 
phase  above  the  surface  of  the  film  to  deposit  desired  materials  onto  a  susceptor. 
CVD  is  a  key  element  in  a  wide  variety  of  industrial  applications,  ranging  from 
the  fabrication  of  microelectronic  circuits,  solar  cells,  and  optical  devices  to  the 
deposition  of  wear  resistant  coatings  onto  high  performance  machine  tools.  In 
a  typical  CVD  reactor,  a  mixture  of  reactants  and  carrier  gas  is  forced  to  flow 
across  a  heated  susceptor.  The  temperature  field  from  the  heated  susceptor 
induces  gas  phase  reactions  to  produce  activated  species  which  then  diffuse  to 
the  surface  reaction  layer  and  decompose  to  produce  a  thin  film. 

The  deposited  films,  whose  thicknesses  range  from  a  few  nanometers  to  a 
few  microns,  must  be  produced  with  controllable  properties  such  as  purity,  com¬ 
position,  thickness,  microstructure,  and  surface  morphology  (see  e.g.,  Jensen  et 
al.  (1991)).  The  tolerance  limits  on  the  properties  of  the  films  vary  with  the 
application.  Currently,  however |  the  required  properties  of  the  majority  of 
industrially  important  thin  films  that  are  produced  by  CVD  have  become  in¬ 
creasingly  difficult  to  achieve.  For  example,  some  most  important  applications 
which  pertain  to  infrared  laser  sources  require  micron  thick  films.  This  man¬ 
dates  epitaxial  deposition  at  rates  higher  than  those  which  can  be  achieved  at 
low  densities.  The  associated  density  gradients  (due  to  large  temperature  gra¬ 
dients  between  the  inlet  and  the  susceptor)  in  a  gravitational  field  will  induce 
natural  convection  flows.  Convection,  in  turn,  influences  the  growth  processes 
in  two  different  ways,  only  one  of  which  is  beneficial.  Convection  increases  mix¬ 
ing  and  the  overall  transport  and,  thus,  the  growth  rate,  which  is  desirable.  On 
the  other  hand,  it  can  also  affect  the  morphology  of  the  solid  adversely.  The 
latter  phenomenon  is  due,  in  part,  to  the  increased  residential  time  resulting 
from  the  slow  diffusional  exchange  of  the  reactant  species  between  the  main 
flow  stream  and  laminar  recirculation  cells.  Thus,  for  both  the  design  of  CVD 
reactors  (including  determining  optimal  operative  conditions  such  as  input  flow 
rates)  and  the  improvement  of  device  properties,  a  qualitative  and/or  quan¬ 
titative  understanding  of  the  transport  processes  in  CVD  reactors  is  of  great 
importance. 

In  the  past  two  decades  steady  progress  has  been  made  in  the  modeling 
of  the  transport  processes  in  CVD  reactors.  For  example,  Moffat  and  Jensen 
(1986),  (1988)  used  a  fully  parabolic  flow  approximation  in  axial  direction  for 
a  three-dimensional  (3D)  horizontal  reactor.  Gokoglu  et  al.  (1989)  studied  the 
deposition  of  Si  with  a  more  detailed  3D  flow  treatment  in  a  similar  reactor 
geometry.  In  Ouazzani  et  al.  (1988),  (1990),  various  2D  and  3D  models  of  a 
horizontal  chemical  vapor  deposition  reactor  were  investigated  and  the  results 
were  compared  with  experimental  data.  Studies  also  have  been  carried  out  to 
investigate  crystal  growth  under  reduced  gravity  conditions  (see  e.g.,  Ouazzani 
et  al.  (1988)  and  Scroggs  et  al.  (1995)).  One  potential  advantage  of  such  con¬ 
ditions  is  that,  under  low  gravity  environment,  buoyancy-driven  convection  is 
reduced  (Ostrach  (1982)).  For  an  in-depth  review  of  the  work  in  these  areas  see 
the  articles  by  Jensen  (1989),  Fotiadis  (1990),  and  Jensen  et  al.  (1991).  Because 
of  the  complexity  in  CVD  reactor  models  (represented  by  systems  of  nonlinear 


Use  of  the  POD  method  in  CVD  reactor 


3 


partial  differential  equations),  the  majority  of  the  models  are  solved  numeri¬ 
cally  by  either  finite- difference  (Coltrin  et  al.  (1984)),  finite- volume  (Ouazzani 
et  al.  (1990)),  or  finite-element  methods  (Jensen  et  al.  (1991)).  On  the  other 
hand,  analytical  investigations  have  used  similarity  transformations  (Pollard 
and  Newman  (1980))  or  separation  of  variable  techniques  (Fujii  et  al.  (1972)). 
These  approaches  are  restricted  to  one- dimensional,  linear,  and  constant  coeffi¬ 
cient  equations.  Consequently,  these  models  neglect  some  of  the  more  important 
nonlinear  effects  in  CVD  processes  such  as  buoyancy,  temperature  dependent 
flow  parameters,  and  nonlinear  coupling  between  the  thermal,  flow,  and  species 
fields.  An  asymptotic  analysis  of  a  CVD  system  which  included  temperature 
dependent  coefficients,  nonlinear  coupling  of  transport  processes  and  Soret  dif¬ 
fusion  effects  was  given  by  Young  et  al.  (1992). 

In  addition  to  the  above  parameter  studies  (effects  of  operating  conditions, 
reactor  geometry,  and  heat  transfer  characteristics  on  flow  patterns  and  growth 
rate  uniformity),  a  more  rigorous  approach  to  the  optimal  design  and  control 
has  also  been  investigated.  In  Ito,  et  al.  (1994),  a  shape  optimization  prob¬ 
lem  with  respect  to  the  geometry  of  the  reactor  and  a  boundary  temperature 
control  problem  were  formulated.  The  material  and  shape  derivatives  of  solu¬ 
tions  to  the  Boussinesq  approximation  were  derived.  Optimality  conditions  and 
a  numerical  optimization  method  based  on  the  augmented  Lagrangian  method 
were  developed  for  boundary  control  of  Boussinesq  flow.  Numerical  calculations 
in  Ito  et  al.  (1995)  indicated  the  effectiveness  of  temperature  control  through 
a  portion  of  the  boundary  for  improving  the  vertical  transport  of  flow  in  the 
cavity. 

Numerical  simulation  has  been  shown  to  be  an  effective  tool  for  the  un¬ 
derstanding  and  improvement  of  CVD  processes.  Of  particular  interest  to  the 
present  investigation  is  the  development  of  a  framework  to  extend  these  simu¬ 
lation  capabilities  into  the  realm  of  optimal  design  and  optimal  control  of  CVD 
reactors.  For  the  optimal  control  of  CVD  processes,  we  begin  by  defining  a  set  of 
process  parameters  which  control  the  system  (e.g.,  flow  rates,  species  concentra¬ 
tion),  and  cost  and  constraint  functions  which  quantify  the  desirable  quantities 
of  the  response  (e.g.,  growth  rate  and  uniform  fluxes  of  reactants  at  the  suscep¬ 
tor).  We  then  perform  the  simulation  of  the  system  process,  and  subsequently 
evaluate  the  cost  and  the  constraint  functions.  These  data  are  then  supplied 
to  some  numerical  optimization  routine  which  modifies  the  process  parameters 
to  reduce  the  cost  functional  value  and  to  satisfy  the  constraints  (Banks,  et 
al.  (1997)).  The  complex  fluid  dynamics  in  CVD  processes  are  described  by 
a  system  of  nonlinear  partial  differential  equations  representing  the  continuity, 
momentum,  energy,  species  and  equation  of  state.  Therefore,  numerical  sim¬ 
ulations  of  such  systems  using  finite  element,  finite  volume,  finite  difference, 
or  spectral  methods  will  lead,  in  general,  to  a  very  large  system  of  ordinary 
differential  equations  rendering  it  inapplicable  in  real  time  estimation  and  con¬ 
trol.  One  approach  to  overcome  these  difficulties  is  to  perform  model-predictive 
control  of  these  distributed  parameter  systems.  Although  this  method  has  the 
merit  of  small  degrees  of  freedom,  they  do  not  represent  the  physical  model,  but 
rather  an  empirical  model  that  is  based  on  input  and  output  of  a  given  system 
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and  may  become  unstable  as  the  operating  condition  of  the  system  changes. 
For  an  efficient  real-time  control  of  CVD  reactors  where  continuous  computa¬ 
tions  must  be  performed,  not  only  the  degrees  of  freedom  of  the  dynamic  model 
must  be  small  but  the  dynamic  model  must  also  be  robust.  In  this  work,  we 
demonstrate  the  feasibility  and  efficiency  of  the  proper  orthogonal  decompo¬ 
sition  technique  (or  the  Karhunen-Loeve  procedure)  in  flow  calculations  and 
the  optimal  control  of  CVD  processes.  The  proper  orthogonal  decomposition 
(POD)  is  a  reduction  method  that  has  been  shown  to  be  an  effective  tool  for 
the  analysis  of  complex  systems  such  as  turbulence  flows,  shear  flows,  pattern 
recognition,  and  weather  prediction  (see  e.g.,  Berkooz,  et  al.  (1993)  and  the 
references  therein).  In  general,  the  discretization  of  nonlinear  partial  differen¬ 
tial  equations  using  finite  element,  finite  volume,  or  finite  different  methods 
involves  basis  functions  that  have  little  to  do  with  the  differential  equation. 
For  example,  piecewise  polynomials  are  used  in  the  finite  element  method,  grid 
functions  are  used  in  the  finite  difference  method,  and  Legendre  or  Chebyshev 
polynomials  are  used  in  some  spectral  methods.  POD,  on  the  other  hand,  uses 
basis  functions  that  span  a  data  set,  collected  from  an  experiment  or  numerical 
simulation  of  a  dynamical  system,  in  a  certain  “optimal”  fashion.  Because  POD 
basis  elements  are  optimal  in  the  sense  that  they  are  the  extractions  of  char¬ 
acteristic  features  of  the  data  set,  only  a  small  number  of  POD  basis  functions 
are  needed  to  describe  the  solution. 

The  organization  of  this  paper  is  as  follow.  We  describe  in  §2  a  horizontal 
CVD  reactor  and  its  mathematical  description.  The  POD  technique  and  its 
mathematical  properties  are  presented  in  §3.  In  §4  we  discuss  issues  concerning 
the  implementation  and  numerical  calculations  using  POD  in  the  context  of 
simulating  the  flow  dynamics  in  the  CVD  process.  Finally,  in  §5  the  method 
is  applied  to  solve  an  optimal  control  problem  to  achieve  film  uniformity  and 
maximum  growth  rate. 

2  Formulation  of  the  Problem 

The  particular  geometry  of  the  organo-metallic  chemical  vapor  deposition  (0M- 
CVD)  teactor  under  consideration  here  features  horizontal  flow  of  the  process 
gases  and  source  vapor/carrier  gas  mixtures  into  an  expansion  section  leading 
into  a  rectangular  channel  that  contains  the  substrate  (see  Figure  1).  The  sub¬ 
strate  wafer  is  mounted  on  a  rotating  induction  heated  SiC  coated  graphite  sus¬ 
ceptor.  The  exhaust  gases  are  vented  through  a  vertical  exhaust  tube.  Loading 
and  unloading  of  substrate  wafers  is  accomplished  through  a  load-lock  chamber 
beneath  the  radio  frequency  (rf)  section  of  the  reactor  that  can  be  evacuated  by 
a  turbomolecular  pump.  After  purging  with  ultra-pure  nitrogen,  sample  trans¬ 
fer  can  be  executed  using  a  magnetic  transfer  rod.  Gas  is  purging  through  the 
gap  between  the  susceptor  and  the  reactor ’s  base  to  avoid  flow  of  gas  mixtures 
to  the  mechanical  workings  behind  the  susceptor.  The  quartz  glass  reactor  is 
connected  at  the  inlet  to  a  source  vapor/process  gas  flow  control  and  switching 
panel  that  directs  individual  streams  of  source  vapor  saturated  carrier  gas  ei¬ 
ther  to  a  vent  line  or  to  the  reactor.  Thus,  pulsed  operation  separating  plugs  of 
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to  load-lock  chamber 


Figure  1:  Schematic  representation  of  a  horizontal,  quartz  reactor  in  a  steel 
confinement  shell 


source  vapor  saturated  carrier  gas  by  plugs  of  high  purity  carrier  gas,  flow  rate 
modulated  flow  or  continuous  flow  can  be  implemented  for  all  source  vapors 
without  change  in  reactor  pressure  or  total  flow.  Two  optical  windows  at  the 
Brewster  angle  of  the  substrate  are  attached  to  the  sides  of  the  reactor.  They 
allow  for  the  real-time  process  monitoring  utilizing  p-polarized  reflectance  spec¬ 
troscopy  (PRS)  (see  e.g.,  Bachmann,  et  al  (1998)  and  the  references  therein). 
This  OMCVD  reactor  has  been  built  in  the  laboratory  of  Prof.  Klaus  Bachmann 
at  North  Carolina  State  University  and  is  now  undergoing  initial  testing. 

To  demonstrate  that  the  POD  method  can  be  implemented  to  approximate 
transport  processes  efficiently  in  the  CVD  reactor,  we  will  restrict  our  study  to 
a  two-dimensional  horizontal  reactor  as  shown  in  Figure  2  (this  can  be  thought 
of  as  a  vertical  “slice”  of  the  actual  3-D  reactor  tube).  Our  study  involving  the 
full  3-dimensional  geometry  will  be  presented  in  subsequent  papers.  We  will 
consider  the  deposition  of  GaN  using  pulsed  trimethyl-gallium  (TMGa)  and 
ammonia  (NH3)  as  source  vapors  and  nitrogen  as  carrier  gas  (see  Figure  3). 
The  function  F(t)  in  Figure  3  represents  the  pulses  of  reactive  gases  entering 
the  reactor.  In  particular,  at  first  only  carrier  gas  flows  through  the  reactor. 
After  the  flow  reaches  steady  state,  a  pulse  of  reactant  (e.g.,  TMGa)  diluted  with 
carrier  gas  enters  the  reactor.  After  the  pulse,  the  reactor  is  then  flushed  with 
carrier  gas.  This  process  is  then  repeated  for  another  reactant.  Furthermore, 
the  following  assumptions  are  made  for  the  mathematical  formulation  of  this 


Inlet 
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process. 


Susceptor 


Figure  2:  Two-dimensional  horizontal  reactor 


AF(t) 


t  t  t  1 

^  Reactant  #  2  =  NH3  Diluted  with  Carrier  Gas 
—  Reactant  #  1  =  TMGa  Diluted  with  Carrier  Gas 

-  Carrier  Gas  =  Nitrogen 

Figure  3:  Chemicals  introduced  into  the  reactor  at  the  inlet 


(i)  Carrier  gas  flows  through  the  reactor  at  all  time. 

(ii)  All  thermo-physical  properties,  conductivity,  viscosity,  mass  diffusivity, 
and  volume  expansion  are  temperature  dependent. 

(iii)  Only  a  trace  amount  of  reactants  mixed  with  carrier  gas  is  allowed  to 
enter  the  reactor  at  each  pulse  so  that  the  steady  state  gas  flow  condition 
and  all  temperature  dependent  parameters  in  (ii)  remain  unchanged. 

(iv)  No  chemical  reaction  takes  place  in  the  gas  phase. 

(v)  Reactions  taking  place  on  the  heated  substrate  are  very  fast  and  equilib¬ 
rium  is  attained  quickly.  Thus,  the  rate  of  GaN  deposition  is  limited  by 
mass  transport. 

(vi)  The  wall  of  the  reactor  is  water  cooled. 

Under  the  above  assumptions,  CVD  processes  can  be  classified  as  a  quasi¬ 
transient  flow  (steady-state  flow  with  transient  species)  and  be  described  by  the 
following  governing  partial  differential  equations  written  in  conservation  form 
(see  e.g.,  Oran  and  Boris  [1987]  and  Oswatitsch  [1956]). 


Outlet 
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Mass  conservation  equation: 


v  •  (pu)  -  0,  (2.1) 

where  p  is  the  density  of  the  carrier  gas  and  u  =  (ft  i ,  U2)  is  the  velocity  vector. 
Momentum  conservation  equation: 

pu  ■  Vu  =  -VP  +  V  •  <7  +  (p  -  po)gj,  (2.2) 

where  the  viscous  stress  tensor  a  has  the  form 

a  —  [/iT(Vu  +  (Vu)2")]  +  (k  -  |jUt)(V  •  u)  •  I. 

Here  we  treat  the  gas  mixture  as  a  Newtonian  fluid  and  P  is  the  pressure, 
Pt  is  the  carrier  gas  viscosity,  k  is  the  bulk  viscosity  which  can  be  neglected 
for  dense  gases  or  liquids,  g  is  the  gravitational  force,  and  po  is  the  reference 
density.  The  ( p  —  po)gj  term  accounts  for  the  natural  convection  effect  caused 
by  the  gravitational  force. 

Energy  conservation  equation: 

cpPu-  VT=  V-(AtVT),  (2.3) 

where  T  is  the  temperature,  and  cp  is  the  heat  capacity  and  A t  is  thermal 
conductivity  of  the  carrier  gas. 

Mass  conservation  of  species: 

dc  1 

-Pu-  Vc=  -V-(^tVc),  (2.4) 

where  c  is  the  mass  fraction  of  TMGa  and  D t  is  the  diffusion  coefficient  of 
TMGa  in  the  carrier  gas.  Here,  without  loss  of  generality,  only  the  transport 
of  TMGa  is  modeled. 

Expansion  of  the  gas  as  it  approaches  the  heated  susceptor  plays  a  major  role 
in  the  flow  behavior  and  it  is  accounted  for  by  using  the  following  Boussinesq 
approximation  for  the  density  as  a  function  of  the  temperature: 


p  =  p0[l-pT(T-T0)],  (2.5) 

where  (3t  is  the  volume  expansion  and  To  is  a  reference  temperature. 

The  boundary  conditions  for  the  above  set  of  equations  (2.2)-(2.5)  are  sum¬ 
marized  in  Figure  4,  where  the  parabolic  profile  for  the  velocity  field  at  the  inlet 
is  described  by 


ui (t)  =  —^2 —y(H  y)  and  u2(t)  =  0, 

with  hmax  =  0.25  m/s  and  H  =  0.535  inch  is  the  height  of  the  reactor.  The 
initial  condition  for  the  species  equation  is  zero  throughout  the  domain  of  the 
rectangular  reactor. 
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Boundary 

conditions 

Velocity 

Temperature 

Species 

Inlet 

Parabolic  profile 

298  K 

F(t) 

Walls 

0 

298  K 

dc/dn  =  0 

Susceptor 

0 

1200  K 

0 

Outlet 

duj  /dn  =  0,  i=l,  2 

dT/dn  =  0 

dc/dn  =  0 

Figure  4:  Boundary  Conditions  for  the  2D  Horizontal  Reactor 


In  this  model,  one  proceeds  first  by  solving  for  the  steady  state  solutions 
of  the  flow  and  energy  equations.  These  solutions  are  then  substituted  into 
equation  (2.4)  to  solve  for  the  time  dependent  species  solution.  Even  though 
the  flow  and  energy  solution  are  decoupled  from  the  mass  transport  analysis,  the 
dynamical  model  is  still  an  infinite  dimensional  system  of  equations.  Standard 
techniques  such  as  finite  difference  or  finite  element  methods  can  be  employed 
to  reduce  the  infinite  dimensional  models  to  finite  dimensional  ones,  but  the 
resulting  degrees  of  freedom  are,  in  general,  too  large  for  practical  considerations 
in  estimation  or  control. 

It  will  be  demonstrated  in  the  sequel  that  a  Galerkin  procedure  employing 
basis  functions  which  are  computed  from  the  proper  orthogonal  decomposition 
can  efficiently  reduce  distributed  parameter  systems  to  low  order  finite  dimen¬ 
sional  dynamical  models  while  maintaining  high  fidelity.  Thus,  this  approach  is 
particularly  suitable  to  treat  optimal  control  and  parameter  estimation  prob¬ 
lems  of  systems  governed  by  partial  differential  equations. 

3  Proper  Orthogonal  Decomposition 

The  proper  orthogonal  decomposition  (POD)  method  has  received  much  atten¬ 
tion  in  recent  years  as  a  tool  to  analyze  complex  physical  systems.  In  principle, 
the  idea  is  to  use  a  reliable  solver  to  produce  a  priori  a  number  of  solutions 
to  the  physical  model  (called  snapshots).  The  POD  technique  is  then  used  to 
produce  an  “optimal”  representation  of  these  snapshots  in  an  “average”  sense. 
Both  notions,  “optimal”  and  “average”,  will  be  made  clear  in  subsequent  dis¬ 
cussions.  The  power  of  the  POD  method  lies  in  its  mathematical  properties 
which  suggest  that  it  is  the  preferred  method  to  use  in  many  applications. 

Proper  orthogonal  decomposition  was  independently  proposed  by  several 
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scientists  including  Karhunen  (1946),  Loeve  (1945),  Pougachev  (1953),  and 
Obukhov  (1954)  (for  recent  surveys  in  this  area  see  the  works  of  Lumley  (1970) 
and  Berkooz  (1992)).  Some  mathematical  theories  of  POD  can  be  found  in 
recent  articles  by  Aubry,  Lian  and  Titi  (1993)  and  Graham  and  Kevrekidis 
(1996).  The  POD  technique  has  been  applied  to  numerous  applications.  One 
such  important  application  was  the  attraction  of  spatial  organized  motions  in 
fluid  flows.  Theodorsen  (1952)  and  later  Townsend  (1956)  observed  and  in¬ 
dicated  that  there  are  large-scale  organized  motions  embedded  in  turbulent 
shear  flows.  Lumley  (1967),  Aubry,  Holmes,  Lumley  and  Stone  (1988),  Sirovich 
(1991),  Berkooz,  Holmes  and  Lumley  (1993),  Berkooz,  Holmes,  Lumley  and 
Mattingly  (1997)  have  adapted  the  POD  technique  to  study  turbulent  flows. 
Other  applications  of  POD  include  channel  flows  by  Moin  and  Moser  (1989), 
Ball,  Sirovich  and  Keefe  (1991),  square-duct  flows  by  Reichert,  Hat  ay,  Biringer 
and  Husser  (1994),  and  shear  flows  by  Rajaee,  Karlson  and  Sirovich  (1993), 
Kirby,  Boris  and  Sirovich  (1990).  Other  scientists  have  also  applied  the  POD 
technique  to  fluid  related  problems.  For  instance,  it  has  been  applied  to  the 
Burgers ’  equation  by  Chambers,  Adrian,  Moin,  Stewart  and  Sung  (1988),  the 
Ginzburg-Landau  equation  and  the  Benard  convection  by  Sirovich  (1989).  Ly 
and  Tran  (1998)  have  used  POD  to  simulate  and  solve  an  optimal  control  prob¬ 
lem  for  Rayleigh-Benard  convection.  Other  interesting  non-fluid  applications  of 
POD  techniques  are  the  characterization  of  human  faces  by  Kirby  and  Sirovich 
(1990)  and  image  recognition  by  Hilai  and  Rubinstein  (1994). 

As  will  be  seen  in  the  following  section,  one  reason  that  POD  is  an  attractive 
method  is  that  it  is  a  linear  procedure.  Its  mathematical  theory  is  based  on 
the  spectral  theory  of  compact,  self-adjoint  operators.  However,  it  should  be 
noted  that  POD  makes  no  assumption  on  the  linearity  of  the  problem  to  which 
it  is  applied  and  this  is  an  extremely  positive  feature  of  this  approach  to  model 
reduction. 

3.1  Mathematical  Aspects  of  POD 

Let  {U  i(x)  :  1  <  i  <  N]  x  E  £2}  denote  the  set  of  N  observations  (also  called 
snapshots )  of  some  physical  processes  over  a  domain  £2.  In  the  context  of  CVD 
process,  these  observations  could  be  experimental  measurements  or  numerical 
solutions  of  velocity  fields,  temperatures,  species  etc.  taken  at  different  phys¬ 
ical  parameters  (Reynolds  number,  input  flow  rates  etc.)  or  time  steps.  The 
POD  technique  is  designed  to  extract  from  this  set  of  observations  a  coherent 
structure,  which  has  the  largest  mean  square  projection  on  the  observations. 
In  other  words,  we  look  for  a  function  4>,  or  the  so-called  POD  basis  element, 
that  most  resembles  {U^(a?))-V1  in  the  sense  that  it  maximizes 

1  N 

-^|(Ui,$)|2,  (3.1) 

i  —  1 


subject  to 


(*,*)  =  \m2  =  i, 
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where  (*,  •)  and  ||  •  ||  denote  the  usual  L2  inner  product  and  L2-norm  over  S2, 
respectively.  We  choose  a  special  class  of  trial  functions  for  to  be  of  the  form: 


N 

$  =  y]aiUil  (3.2) 

i  —  1 

where  the  coefficients  cii  are  to  be  determined  so  that  given  by  the  expression 
(3.2)  provides  a  maximum  for  (3.1).  To  this  end,  let  us  define 

1  *  f 
K(a?,x/):= — and  R<F  :=  /  K  ( x ,  xf )  ( ,U)  dxl , 

N  Jn 


where  R  :  L2(D)  -►  L2(D). 

Then  straightforward  calculations  reveal  that 

(R$,$)  =  /  R$(x)$(x)dx 

Jn 

—  /  /  K (x,  xf)<f>(xf)dxf<f>(x)di 

Jn  Jn 

N 


=  V  /  [  Ui(x)Ui(x,)^(x,)dx,^(x)da 

^  •_!  JnJn 


Furthermore,  it  follows  that 

(R$,^)  =  (<F,R^)  for  any  GL2. 

Thus  R  is  a  nonnegative  symmetric  operator  on  L2(f2).  Consequently,  the  prob¬ 
lem  of  maximizing  the  expression  (3.1)  amounts  to  finding  the  largest  eigenvalue 
to  the  eigenvalue  problem 


R<F  =  A<F  subject  to  ||<F||  =  1,  (3-3) 


or 


/ 

Jn 


K(a?,  xr)Q{xr)dxf  —  A<F  with  ||4>||  =  1. 


(3.4) 


Substituting  expression  (3.2)  and  the  definition  of  K  into  equation  (3.4),  we 
obtain 


N  N 


=  1  k- 1 


If  N 

(—  /  \Ji{x!)\Jk{xf)dxf)ak]Ui(x)  =  XajUi(x). 

^  Jn 


This  can  be  rewritten  as  the  eigenvalue  problem 


CV  =  AV 
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where 


Oik  — 


1  f 

=  —  I  XJi(x)XJ k(x)dx  and 

N  Jn 


V  = 


ai 

aN 


(3.5) 


Since  C  is  a  nonnegative  Hermitian  matrix,  it  has  a  complete  set  of  orthogonal 
eigenvectors 


V1  = 


with  the  corresponding  eigenvalues  Ai>A2>***>Aj\r>0.  Thus,  the  solution 
to  the  optimization  problem  for  (3.1)  is  given  by 

N 


« 1 

«1 

Ki 

«2 

< 

to 

II 

<*2 

< 

II 

a? 

.  aJV  . 

_  °N  _ 

- 1 

.  ^ 
G 

_ i 

*1  -  E  a]  U, 


i  —  1 


where  a)  are  the  elements  of  the  eigenvector  V1  corresponding  to  the  largest 
eigenvalue  Ai.  The  remaining  POD  basis  elements,  <^,  i  =  2,  .  .  .  ,iV,  are  ob¬ 
tained  by  using  the  elements  of  other  eigenvectors,  V*,  i  =  2,  .  .  . ,  N.  Moreover 
using  the  orthogonality  of  {V^  :  1  <  k  <  N }  and  the  imposed  condition: 


N 


v*  .  v*'  =  ^2 


a)  a)  —  i  N  A*./ 


8  =  1 


k  =  k' 
k  ±  k', 


we  obtain 


=  /  $>k(x)$k>(x)dx 

Jn 

.  n  i v 

-  /  Ea»'U<^)Eaf  \Jj(x)di 

Jn  8  =  1  j  =  1 

IV  IV 

=  E«‘A'E<v  yu«(*)u, ■<*)<**)“*' 


*= i  j=i 

JV  JV 

E«8^Ec 

*=i  j=i 

JVV*  •  cv*' 

iVV*  •  Ajfe/V*' 
NXkfVk  •  V*' 
f  1  Jb  =  A' 

\  0  Ar  ^ 


7ft' 
»J  i 
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Thus,  the  POD  basis  {$1,  $2,  .  .  . ,  forms  an  orthonormal  set. 

Remarks.  An  alternative  approach  for  finding  the  solution  to  maximization 
of  (3.1)  is  by  using  the  so-called  Rayleigh- Ritz  method  for  finding  eigenvalues 
(see  Kanwal  (1997)  p.  176).  Since 

N  N  N  N 

(R<£,  =  (Rfr:  akD$k)  =  R  ik&i&k  ? 

i-l  k- 1  i-l  k- 1 

and 

N  N 

\m2  =  (*,*)  =  Y,Y,NCikaiak> 

i-l  k- 1 

where  Rik  =  (RU^U*)  =  (Uj,RUfc)  and  Cik  —  ^(Uj,Uj,)  =  -pL-V  V)  as 
above  in  (3.5),  the  problem  of  maximizing  (3.1)  can  be  transformed  into  an 
extremal  problem  in  multi- variable  calculus  with  {ai,  02,  ,  ajv}  as  variables. 

Recalling  the  method  of  Lagrange  multipliers,  we  define 

N  N 

G(ai,  .  .  ,,aN)  :=  ^  ^  Rikaiak  -  XNCikaiak, 
i  =  l  k= 1 

where  A  is  the  Lagrange  multiplier.  Equating  $G/'&&{  to  0  (a  necessary  condi¬ 
tion  for  optimality)  we  find 

N 

^(Rji  -  XNCik)ak  —  0  for  i  =  l,...,N.  (3.6) 

k  =  l 

Equation  (3.6)  has  a  nontrivial  solution  if  and  only  if  the  determinant  of 
(R  —  AAC)  =  0.  This  determinant,  when  expanded,  yields  an  Ath  degree  poly¬ 
nomial  in  A  which  can  be  shown  to  have  A  nonnegative  real  roots,  but  not  neces¬ 
sarily  distinct  (these  can  be  ordered  in  descending  order  as  Ai  >  A2  >  •  •  •  >  Ajy) 
(Kanwal  (1997)).  Moreover  by  multiplying  equation  (3.6)  by  ai  and  summing 
over  i  =  1,  .  .  . ,  A,  we  obtain  A  =  (R<1>,  <£).  It  can  be  shown  that  solution  pair 
(A^,  V%)  where  V*  =  (a\ ,  .  .  . ,  a^),  is  the  eigenvalue  and  eigenvector  of  the  ma¬ 
trix  C.  One  then  uses  the  definition  of  $  in  (3.2)  to  obtain  {$1,  $2,  .  .  . ,  ^jv} 
and  uses  the  fact  that  the  matrix  C  =  [C^]  is  Hermitian  to  justify  the  or- 
thonormality  of  {<l>£}’s  as  in  previous  argument. 

3.2  Optimality  of  the  POD  Basis 

Suppose  that  we  have  a  signal  v(x,t)  with  v  £  L2(S2,  [0,T])  and  an  approxi¬ 
mation  of  vN  of  v  with  respect  to  an  arbitrary  orthonormal  basis  ipi(x),  i  = 
1,2,...,  A: 

N 

vN(x,t)  —  ymqm*)- 

i  —  1 
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If  the  ipi(x)  have  been  nondimension  alized,  then  the  coefficients  cti  carry  the 
dimension  of  the  quantity  vN .  If  vN(x,  t)  denotes  the  velocity  and  <  •  >  denotes 
the  time  average  operator,  then  the  average  kinetic  energy  per  unit  mass  is  given 

by 


f  N 

<  /  vN(x,  t)vN*(x,t)  dx  >=<  aj{t)at  ( t )  >  . 

,= i 

Consequently,  the  expression  <  >  represents  the  average  kinetic  energy  in 

the  ifb_mode.  The  following  lemma  establishes  the  notion  of  optimality  of  the 
POD  method. 

Lemma  3.1  Let  {$i,  $2,  •  •  • ,  $jv}  denote  the  orthonormal  set  of  POD  basis 
elements  and  (A] ,  A2,  .  .  . ,  AjyJ  denote  the  corresponding  set  of  eigenvalues .  If 

N 

vN(x,t)  - 

i  —  1 

denotes  the  approximation  to  v  with  respect  to  this  basis ,  then  the  following 
hold: 

(a)  <  bi(t)bj(t)  >=  6ij \j  (that  is,  the  POD  coefficients  are  uncorrelated). 

(b)  For  every  N , 

N  N  N 

Fi  <  bi(t)b*  w  >=  Fi Xi  -  H  <  a*coa*co  >  • 

i — 1  i  —  1  i  —  1 

The  proof  of  this  lemma  is  straight  forward  from  the  optimality  of  the  eigen¬ 
values  and  can  be  found  in  Berkooz,  Holmes,  and  Lumley  (1993).  This  lemma 
establishes  that  among  all  linear  combinations,  the  POD  is  the  most  efficient, 
in  the  sense  that,  for  a  given  number  of  modes,  N ,  the  projection  on  the  sub¬ 
space  used  for  approximation  will  contain  the  most  kinetic  energy  possible  in 
an  average  sense. 

3.3  Model  Reduction  Features  of  POD  Approximations 

To  this  point  we  have  not  discussed  any  model  reduction  features  associated 
with  using  POD  basis  elements  in  approximation  schemes.  In  the  construc¬ 
tion  described  above,  the  number  N  may  be  large,  100  —  1000  or  even  more, 
depending  on  the  complexity  of  the  dynamics  represented  in  the  “snapshots” 
U.?; .  In  general,  one  should  take  N  sufficiently  large  so  that  the  snapshots  Ux- 
contain  all  salient  features  of  the  dynamics  being  investigated.  Thus,  the  POD 
basis  functions  <^,  used  with  the  original  dynamics  in  a  nonlinear  Galerkin 
procedure,  offers  the  possibilities  of  achieving  a  high  fidelity  model,  albeit  with 
perhaps  a  large  dimension  N . 

To  achieve  model  reduction,  one  chooses  M  <C  N  and  carries  out  a  nonlin¬ 
ear  Galerkin  procedure  with  the  set  of  elements  {$1,  $2,  .  .  . ,  The  crucial 
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question  is  how  to  choose  M.  As  we  indicated  in  the  previous  section,  YliLi  ^ 
represents  the  average  kinetic  energy  contained  in  the  first  M  modes  and  hence 
to  capture  most  of  the  energy  of  the  system  contained  in  the  N  POD  ele¬ 
ments,  it  suffices  to  choose  M  so  that  J2iLi  ~  'V:  •  Indeed,  the  ratio 

1  /  J2iLi  yi^ds  the  percentage  of  the  total  kinetic  energy  in  the  N 

POD  elements  that  is  contained  in  the  first  M  POD  elements.  Since  the  asso¬ 
ciated  POD  eigenvalues  are  ordered  Ai  >  A2  >  •  •  •  >  Ajy,  one  can  reasonably 
expect  to  achieve  a  high  percentage  of  the  total  kinetic  energy  in  a  reduced 
model  of  order  M  with  M  sufficiently  smaller  than  N .  For  the  CVD  examples 
detailed  below,  the  POD  system  was  constructed  for  N  =  200  and  a  reduced 
order  model  with  M  =  10  yielded  a  ratio  of  .999,  resulting  in  a  truly  signifi¬ 
cant  computational  savings  over  the  finite  element  model  (2,400  quadrilateral 
elements)  of  dimension  14,801  used  to  generate  the  200  snapshots. 


4  Application  of  POD  for  the  Simulation  of  CVD 

Processes 


We  return  to  the  CVD  example  of  §2  to  apply  the  POD  techniques.  Under 
assumption  (iii)  of  §2,  that  is  only  a  trace  amount  of  TMGa  mixed  with  carrier 
gas  is  used,  the  steady  state  flow  and  energy  solutions  can  be  decoupled  from 
the  mass  transport  equation.  Therefore,  we  first  solve  equations  (2.1)-(2.3) 
and  (2.5)  for  the  steady  state  solutions  of  velocity  field,  temperature,  density 
and  pressure.  These  solutions,  depicted  graphically  in  Figure  5,  are  computed 
using  a  commercial  fluid  dynamics  package  called  FIDAP  version  7.6  which 
employs  the  finite  element  method.  In  our  simulations  we  used  2,400  (9-nodal 
quadratic)  quadrilateral  elements.  Here  and  in  all  subsequent  plots,  higher 
numerical  values  are  represented  by  darker  shaded  region. 

Once  the  flow  has  reached  steady  state  condition  in  the  reactor,  we  introduce 
TMGa  into  the  reactor  from  the  inlet.  As  mentioned  earlier,  the  term  F(t) 
used  in  the  boundary  condition  for  the  species  at  the  inlet,  c(<r,  2)|*nje|.  =  T(tf), 
describes  the  incoming  pulses  of  TMGa.  Using  the  steady  state  velocity  and 
density  solutions  (for  the  terms  u  and  p)  and  the  steady  state  temperature 
solution  (to  compute  temperature  dependent  mass  diffusivity  Dt )  in  equation 
(2.4),  we  simulated  the  species  equation  for  zero  initial  condition  and  various 
boundary  conditions  F(t)  using  FIDAP.  Figure  6  displays  TMGa  profiles  at 
different  time  steps  from  0.05  to  0.75  seconds  corresponding  to  the  boundary 
condition  c(x,  t)|*nje|-  =  F(t)  =  1,  for  all  t  >  0.  We  note  that,  depending  on  the 
duration  of  TMGa  introduced  at  the  inlet,  different  TMGa  profiles  will  result  in 
the  reactor.  In  particular,  for  the  following  boundary  condition  (see  also  Figure 
8): 


^4tnlet  =  {  0^  1 


0  <  t  <  I  in. 

Tin  <  *  -  Tin+Tout . 


Figure  7  depicts  two  different  reactant  profiles  at  different  time  steps  from  0.05 
second  to  Tou^  =0.4  second  corresponding  to  Tin  =  0.3  second  and  Tin  = 
0.75  second.  In  the  sequel,  we  will  demonstrate  that  the  proper  orthogonal 
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Steady-State  Flow  with  Nitrogen  as  Carrier  Gas 


CO 

c 

Cl) 

Q 


Figure  5:  Steady  state  velocity  field,  temperature,  density  and  pressure 


decomposition  method  implemented  in  a  Galerkin  procedure  can  be  used  to 
simulate  both  accurately  and  efficiently  the  complicated  species  dynamics  for 
different  boundary  conditions  (i.e.,  different  Tjn  values). 

4.1  Construction  of  POD  Basis  Vectors 

Consider  the  following  semi-discrete  nonlinear  equation  of  the  form 

^  =  g(t,w(t)),  for  ten,  wex,  g-.RxX^X,  (4.1) 

where  X  is  a  finite-dimensional  space.  If  finite  element  procedures  were  used  to 
obtain  this  semi-discrete  problem,  then  the  choice  for  X  would  be  defined  by 
spanj^i,  (f2 ,  .  .  . ,  where,  for  example,  tpi  are  piecewise  polynomial  functions 
(e.g.,  the  so-called  hat  functions,  or  spline  functions).  In  the  POD  technique, 
however,  we  will  make  a  different  choice  for  the  approximating  space.  Let  XPOD 
denote  the  POD  space  such  that  XPOD  C  X.  The  procedure  for  computing 
XPOD  consists  of  the  following  steps. 
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Species  Introduced  into  Reactor 


d  i 


Figure  6:  TMGa  profile  due  to  the  continuous  input  of  reactant  at  the  inlet 


(i)  Obtain  the  snapshots.  We  first  allow  the  chemical  mixture  of  TMGa  and 
N2  to  enter  the  reactor  at  the  inlet  for  0.5  second.  After  0.5  second,  we 
shut  off  TMGa  and  allow  only  carrier  gas,  Jf|,  to  flow  in  for  one  sec¬ 
ond,  during  which  we  solve  the  species  equation  at  200  time  steps  (then 
snapshots)  {ci(x),  C2(x),  .  .  . ,  C2oo(x)}  at  an  increment  of  0.005  second  for 
x  E  £2  (here,  denotes  the  two-dimensional  rectangular  domain  as  de¬ 
picted  in  Figure  2).  These  snapshots  are  pointwise  discrete  data  of  species 
over  f2,  which  have  been  computed  using  FIDAP.  Some  sample  snapshots 
of  species  are  displayed  in  Figure  9.  We  note  that  at  time  t  —  0.8  sec¬ 
ond  (see  snapshot  #  160)  most  of  the  TMGa  has  been  carried  out  of  the 
reactor. 
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Figure  7:  TMGa  profiles  corresponding  to  Tjn  equals  to  .3  second  (left)  and 
0.75  second  (right) 


(ii)  Compute  the  covariant  matrix  C.  The  matrix  elements  of  C  are  given  by 

c«  =  2M  J^C{(*)Ck&)dx' 

for  z,  k  =  1,2,...,  200. 

(iii)  Solve  the  eigenvalue  problem  CV  =  AV.  We  recall  that  since  C  is  a 
nonnegative,  Hermitian  matrix,  it  has  a  complete  set  of  orthogonal  eigen¬ 
vectors 
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with  the  corresponding  eigenvalues  arranged  in  ascending  order  as  Ai  > 
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Figure  8:  Schematic  diagram  for  the  boundary  condition  of  TMGa  at  the  inlet 


(iv)  Compute  the  POD  basts  vectors.  The  POD  basis  elements  <F?;(x)  such  that 
XPOD  =  span{<Fi(x),  <F2(x)j .  .  . ,  4>2oo(x)}  are  defined  as 

200 

=  ^  ^  ai  ci $ 
i-1 

where  1  <  k  <  200  and  af  are  the  elements  of  the  eigenvector  corre¬ 
sponding  to  the  eigenvalue  A*. 

4.2  Reconstruction  of  Solutions  Using  POD  Basis  Vectors 

We  next  consider  the  species  equation 
dc  1 

—  +  u  •  Vc  =  -V  •  (pJPt  V c) ,  c(x,  0)  =  flf(x),  (4.2) 

where  u,  and  the  temperature  dependent  parameters  /?,  and  Dt  are  obtained 
from  the  steady  state  solutions  of  the  coupled  system  (2.1)-(2.3)  and  (2.5)  with 
boundary  condition  as  described  in  Figure  4.  In  this  section,  we  will  consider  the 
problem  of  approximating  the  infinite-dimensional  equation  (4.2)  by  a  sequence 
of  finite-dimensional  problems  using  a  combination  of  Galerkin  approximations 
and  POD  basis  elements.  We  first  formulate  the  species  equation  (4.2)  in  a 
variational  or  weak  form;  that  is,  we  seek  a  solution  t  — *  c(t)  on  0  <  t  <  T, 
with  c(t)  E  Ff1(S2)  satisfying 

(ct ,  ip)  +  (u  *  Vc,  ip)  =  (— ^-Vc*  Vp,  ip)  —  (Dt Vc,  Vip)  +  f  DTipVc-nds  (4.3) 
P  Jim 

for  all  'ill  E  Ff1(S2),  along  with  initial  condition 

c(0)  =  g. 

Here,  n  denotes  the  unit  outward  normal  vector  to  the  boundary  dfl  of  the 
domain  Q. 
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Figure  9:  Snapshots  #1,  40,  80,  120,  160  and  200 


As  we  have  discussed  in  the  previous  section,  fewer  than  200  POD  basis 
vectors  are  used  to  approximate  the  species  solution  thereby  achieving  a  reduced 
order  model.  That  is,  the  energy  contained  in  the  first  M  POD  modes,  M  <C 
200,  represents  most  of  the  energy  in  the  system,  if  we  choose  M  so  that 

M  200 

i  —  1  i  —  1 

In  our  case,  we  found  that  the  first  10  POD  basis  functions,  displayed  in  Figure 
10,  capture  over  99.9%  of  the  characteristics  of  the  200  observations.  Here,  for 
i  =  1 .....  1  (J,  POD  #  i  refers  to  the  graph  of  the  function  $*(x)  for  x  E  Q. 
More  precisely, 

A- 

^  -  goo 

—i  A 

Therefore,  it  is  reasonable  to  approximate  the  species  solution  by  using  only 
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The  first  of  these  equalities  results  from  using  (4.4)  in  a  strong  form  of  (4.2), 
multiplying  by  elements  and  integrating,  while  the  equivalent  second  equality 
results  from  using  (4.4)  directly  in  the  weak  form  (4.3)  in  the  usual  Galerkin 
procedure.  The  solutions  to  the  initial  value  problem  (4.5)  yield  the  coefficients 
of  the  POD  basis  function  approximation  (4.4). 


4.3  Simulation  Results 


In  this  section,  we  examine  the  accuracy  and  efficiency  of  the  low-dimensional 
dynamical  model  (4.5)  obtained  from  the  Galerkin  procedure  employing  POD 
basis  functions  by  comparing  its  solution  with  the  solution  obtained  from  the 
fluid  dynamics  package  FIDAP.  More  specifically,  in  FIDAP  the  domain  £2 
is  discretized  using  2400  (9-nodal  quadratic)  quadrilateral  elements.  Conse¬ 
quently,  by  using  FIDAP,  a  system  of  14,801  ordinary  differential  equations  has 
to  be  solved  for  the  coefficients  of  the  basis  function  approximation.  We  point 
out  again  that  in  our  formulation  using  a  Galerkin  procedure  with  POD  basis 
functions,  the  resulting  approximation  is  a  system  of  10  ordinary  differential 
equations  for  the  coefficients  of  the  POD  basis  function  approximation  (4.4). 

Figure  11  compares  the  reduced  solution  using  10  POD  basis  functions  to 
the  full  solution  obtained  from  FIDAP  employing  2400  (9-nodal  quadratic) 
quadrilateral  finite  elements.  Specifically,  both  solutions  are  obtained  with 
boundary  condition 


c(^)linlet  =  {  1 


0  <  t  <  Tin, 

T^  <  ^ 


where  X*n  =  0.5  second.  Qualitatively,  the  reduced  POD  basis  results  agrees 
favorably  with  the  full  FIDAP  calculations.  Quantitatively,  the  upper  plot  in 
Figure  12  graphs  the  L 2  norms  of  FIDAP  solutions,  ||cfidap(*j  ^)|U2(n),  and  the 
L2  absolute  errors  between  POD  solutions  and  FIDAP  solutions  ||cpod(*G)  — 
CFiDAp(*G)|U2(n)-  The  bottom  plot  in  Figure  12  graphs  the  L2  norms  of  FIDAP 
fluxes,  ||^rFiDAp(*G)||L2(n)>  and  the  L2  absolute  errors  between  POD  fluxes  and 
FIDAP  fluxes  ||^7pod('^)  —  ^FiDAp('G)|U2(ns)  Here,  the  flux  of  the  species 
above  the  susceptor  (fis)  that  dictates  the  rate  at  which  species  arrives  at  the 
surface  is  defined  by 


f(Z,t)  :=  (£»TVc(x,t)  +  pc(x,t) u)  •  n|Sens- 

The  flux  of  species  above  the  susceptor  is  used  in  the  next  section  to  define  the 
cost  functional  in  an  optimal  control  problem.  Therefore,  it  is  important  that 
this  quantity  can  be  computed  accurately.  For  a  comparison,  Figure  13  depicts 
the  fluxes  at  the  susceptor  computed  using  POD  basis  functions  and  FIDAP 
software  package.  Both  results  agree  remarkably  well. 


0.525  s  0.42  s  0.315  s 
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After  Species  in  for  0.500  sec  in  the  Reactor 


After  Species  in  for  0.500  sec  at  the  Susceptor 


Figure  12:  Top:  the  L2-norm  of  FIDAP  solution  (o),  ||cfidap(*j^)||,  and  the 
^2-norm  error  (*),  ||cpod(*3  t)  ~  cfidap(-, Oil*  Bottom:  the  L2-norm  of  FIDAP 
flux  (o),  1 1 j^fidap (* ?  1 1 j  and  the  ^2-norm  of  the  flux  error  (*),  ||J7pod(*j^)  — 

^FIDAP  ('5 1)  1 1 
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FIDAP  FLUX  (dotted)  -vs-  POD  FLUX  (dashed)  at  the  Susceptor  *  1 .0e-06 


Figure  13:  Fluxes  of  species  at  the  susceptor  using  FIDAP  (.)  and  POD  basis 
functions  (-)  at  different  times 

We  recall  that  POD  basis  functions  were  computed  from  the  snapshots  cor¬ 
responding  to  the  boundary  condition  (4.3)  with  T*n  =  0.5  second.  To  be  useful 
in  developing  both  open  loop  and  feedback  control  strategies,  we  would  like  for 
these  same  reduced  number  of  POD  elements  to  yield  good  approximations 
under  other  flow  input  conditions.  That  is,  we  would  like  to  demonstrate  that 
using  exactly  these  10  POD  modes,  we  are  able  to  construct  POD  solutions 
for  different  durations  of  species  introduced  into  the  reactor.  Particularly,  we 
solved  the  same  system  (4.5)  with  initial  conditions  ay  (0)  =  (c(x,  Tin),  <F)  for 
Tin  =  0.3  second  and  Tin  =  0.75  second.  Figure  14  depicts  the  L2-norms  of  the 
solution  errors  and  iv2-norms  of  the  flux  errors  above  the  susceptor  using  POD 
basis  functions  and  FIDAP  for  T*n  =  0.30  second  and  T*n  =  0.75  second. 
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After  Species  in  for  0.300  sec  in  the  Reactor 


At  the  Susceptor 


After  Species  in  for  0.750  sec  in  the  Reactor 


Figure  14:  (left  columm,  Y*n  =  0.3  second,  right  columm,  X*n  =  0.75  sec¬ 
ond):  Top,  the  L2-norm  of  FIDAP  solution  (o),  |)cfidap(*>  Oil*  an<^  the  solution 
error  (*),  ||cpoD(*,t)  —  cfidap(-,  t)\\]  Bottom,  the  L2-norm  of  FIDAP  flux  (o), 
||^FiDAp(*,t)||,  and  the  L2-normofthe  flux  error  (*),  ||JFPOD(.,t)- .FfidapO,  t)\\ 


5  An  Optimal  Control  Problem 

In  this  section  we  demonstrate  the  use  of  reduced  POD  models  in  an  open  loop 
optimal  control  problem.  While  this  example  involves  only  one  control  variable, 
it  does  illustrate  effectively  the  potential  of  POD  methods  in  control  problems. 

In  the  case  of  GaN  heteroepitaxy  film  growth  employing  pulsed  trimethyl- 
gallium  (TMGa)  and  ammonia  (NH3)  as  source  vapors,  depending  on  the  delay 
between  the  TMGa  and  NH3  source  vapor  pulses,  carry-over  of  TMGa  frag¬ 
ments  from  one  precursor  pulse  cycle  to  the  next  may  occur.  This,  in  turn, 
establishes  a  surface  reaction  layer  (SRL),  consisting  of  mixture  of  reactants 
and  products  of  the  chemical  reactions  that  drive  the  epitaxial  growth  process. 
The  thickness  and  composition  of  the  SRL  depends  on  the  relative  heights  and 
widths,  i.e.  T*n,  of  the  employed  TMGa  and  NH3  source  vapor  pulses  and  their 
repetition  rate.  More  specifically,  if  T*n  is  small,  then  the  flux  above  the  sus¬ 
ceptor  is  uniform.  Yet  this  will  take  a  long  time  to  grow  a  film  which  makes 
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it  inviable  in  industrial  applications.  On  the  other  hand,  by  allowing  massive 
chemical  input  at  the  inlet  (i.e.,  T*n  is  large)  we  will  speed  up  the  growing 
process;  but  the  flux  above  the  susceptor  may  not  be  uniform.  Therefore  it 
is  desirable  to  perform  an  optimization  procedure  to  determine  the  most  de¬ 
sirable  T*n  to  achieve  both  film  uniformity  and  fastest  possible  growth  rate. 
Mathematically,  we  would  like  to  find  the  optimal  T*n  so  that  the  flux  fluctu¬ 
ation,  ',  is  small  while  at  the  same  time  we  maximize  the  flux,  T ,  to  the 
substrate.  Consequently,  one  way  to  formulate  a  quantitative  representation  of 
these  conflicting  desires  is  to  seek  T*  so  that  the  cost  functional 


a; 

Cost  =  I  -f - , 

J 


i 

-  5.0  x  1012’ 


(5.1) 


is  minimized,  where 


/'/in  +  /out  f  d  l2 
1=1  /  \^-F\  dxdt, 

Jo  J ns  Ox 

and 

/Tin+Tout  f  , 

J  =  /  /  \r\2dxdt. 

tl  0  J  £2  s 

Here,  Tou^  denotes  the  time  delay  between  the  source  vapor  pulses  and,  based  on 
experimental  results,  we  take  Tou^  «  0.600  second.  The  above  cost  functional 
(5.1)  is  minimized  subject  to  the  system  (2.2)- (2.5)  along  with  their  boundary 
conditions.  The  constant  a  in  (5.1)  is  scaled  to  make  I  and  j  proportioned.  It 
also  bears  the  trade-off  value  which  represents  on  one’s  desires  to  balance  film 
uniformity  and  faster  growth  rate. 

Given  a  value  of  T*n  the  system  of  equations  (2.2)-(2.5)  is  solved  for  the 
species  solution  which  is  then  substituted  into  the  formula  (4.3)  for  the  flux, 
T .  Hence,  the  above  optimization  problem  is  an  unconstrained  minimization 
problem  involving  one  parameter,  fT*  .  We  carried  out  computations  for  such  a 
control  example.  For  these,  we  chose  the  DUVMGS  subroutine  in  the  Fortran 
IMSL  library  [1989]  which  is  designed  to  minimize  a  nonsmooth  function  with 
double-precision  accuracy.  The  subroutine  DUVMGS  uses  the  so-called  golden 
section  method  to  search  for  its  minimal  point.  Other  choices  of  minimization 
techniques  are  possible.  We  found  that  with  the  value  of  a  given  in  (5.1),  the 
optimal  width  of  the  pulses,  T-*n,  is  0.395  second.  The  time-dependent  fluxes 
above  the  susceptor  for  the  optimal  are  compared  in  Figure  15  with  non- 
optimal  fluxes  computed  using  T*n  =  0.300  second  and  X*n  =  0.650  second.  We 
note  that  Tjn  must  be  at  least  .300  second  which  is  the  fastest  on-off  switch  in 
the  flow  control  panel.  The  flux  above  the  susceptor  for  T*n  =  0.300  second  is 
uniform,  i.e.,  the  curve  is  flat,  however,  the  growth  rate  is  small.  On  the  other 
hand,  for  X*n  =  0.650  second,  the  growth  rate  speeds  up,  but  non-uniformity 
of  film  growth  results.  The  optimal  solution  lies  between  the  above  two  non- 
optimal  curves.  This  solution  which  corresponds  to  the  optimal  duration  of  the 
source  vapor  pulses  achieves  both  film  uniformity  and  fastest  possible  growth 
rate. 
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Flux  above  the  Susceptor  *1 .0e06 
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TJn  =  0.300  s  T*_in  =  0.395  s  (solid),  TJn  =  0.650  s  ( — ) 


Figure  15:  Comparing  the  optimal  flux  solution  with  two  non-op timal  flux 
solutions 
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